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ABSTRACT 


In this thesis, the Always-Convergent Iterative Noise 
Removal and Deconvolution Method of Ioup is applied as a 
single-filter in the transform domain to deconvolution with 
both narrow and wide Gaussian impulse reponse functions. The 
wraparound error for both cases is also studied. A method is 
developed by which one can find the optimum iteration number 
for single-filter iterative deconvolution of sampled data. 
The method employs the mean square error (MSE), the square 
of ohe difference between the deconvolved result and the 
inpuo, for optimization. The MSE decreases as the 

deconvolution iterations proceed, but at the optimum 
iteration number, the MSE starts to increase. This procedure 
is repeated for signal- to-noise ratio of 10 to 150. The 

optimum iteration number and the MSE are plotted vs SNR. By 
knowing the SNR for a particular experiment, one can find 

the optimum iteration number and MSE. 

The result shows tnat the narrow impulse response has 

smaller optimum iteration numbers and MSE than the wide 

Gaussian, which means the deconvolution of the narrow 

Gaussian has better results than the wide Gaussian. The 
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optimum iteration numbers for the narrow Gaussian range from 
2.8 to 19*2 and for the wide Gaussian from 3*0 to 53«0. 


This thesis also shows that even when noise 
zeros can be added to the original function t' 
wraparound error. An optimum number of zeros 
chosen to save computer time. 


is present, 
reduce the 
has to be 
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CHAPTER I 


THEORETICAL BACKGROUND 


1 . 1 INTRODUCTION 


For much signal processing work, the main task is to recover the 
input signal from -die output signal and the impulse response function. 
Fcr many linear systems, the relation between the output and input may 
be given by a convolution operation. In other words, the convolution 
of The input with the impulse response function of the system gives 
the output. Thus, in order to recover the input signal, one has to 
deconvolve the output signal with the impulse response function. 

There are many deconvolution methods. The iterative deconvolution 
me uiods are among those that are widely used. 

In his thesis, an input signal with three peaks is chosen to be 
convolved with a narrow and a wide Gaussian impulse response function. 
Gaussian type noise is added to the corresponding output. The 
Always-Conv ergent Iterative Deconvolution Method of Ioup (Ioup, 1981) 
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is used "Do deconvolve the output. The mean square error (MSE) between 
the deconvolved result and the input is calculated for the 
optimization. The MSE decreases as the deconvolution iterations 
proceed, but at the optimum iteration number, the MSE starts to 
increase. The procedure is repeated for signal- to-noise ratios of 10 
to 150. By knowing the optimum iteration numbers and mean square 
errors for the narrow and wide Gaussian cases, one is able to evaluate 
the deconvolution for these two cases - 

The deconvolutions are also performed at different data lengths for 
the same signal- to-noise ratios to check the wraparound error . 
Conclusions about whether one need to add zeros in the time domain 
(extending the period) to reduce the wraparound error are made after 
the deconvolutions of different data lengths are done. 

This thesis contains four chapters. Chapter One gives a brief 
mathematical background, including Fourier Transforms, convolution 
and deconvolution, and especially iterative deconvolution methods; 
Chapter Two introduces the Always-Convergent Deconvolution Method of 
Ioup, with discussion and conditions for using the method; Chapter 
Three describes Gaussian type noise and shows how to generate the 
noise; In Chapter Pour, the Always-Convergent Method is applied to 
deconvolve synthetic data, and conclusions about the results and the 
conditions for the deconvolution are made. 
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1 .2 THE PDURIER TRANSFORM 


The mathematical expressions of the Fourier Transform and its 
inverse are: 


F(s) - J f(x) exp(-i2llxs) dx (1.1) 

and 


^( x ) _ J F(s) exp(+i2IIxs) ds, 


(1 . 2 ) 


where F(s) is the forward Fourier Transform of f(x) and f(x) is the 
inverse Fourier Transform of F(s) (Bracewell, 1984). This process 
is possible if f(x) satisfies the following conditions: 

1 ) The integral of from minus infinity 
"Do plus infinity of |f(x)j exists. 

2) Any discon tinuties in f(x) are finite. 

For discrete data, the Fourier Transforms become the following 
summations: 


and 


oc 

F(n) = ^ f (m) exp(-i2Hnm) 

w * ~*c 

f(m) = £. F(n) exp(+i2llnm) , 


(1.3) 

(1.4) 
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inhere m and n are discrete variables. 


1 .3 CONVOLUTION AND THE CONVOLUTION THEOREM 


1*3.1 CONVOLUTION 

For a linear system The input convolved with the impulse response 
function of the system gives the output (Robinson, I960). The 
following diagram shows a convolution model of a linear 
shift- invariant system: 


system 



The mathematical expression for the convolution is 
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(1.5) 


h(x) = f f(u) g(x-u) du = f(x)*g(x) 


where h(x) is the output signal, f(x) is the input signal and g(x) is 
"the impulse response function of the instrument. 

For the discrete situation, the convolution may be given by the 
following equation: 

n n = ^m Sn-m> (1*6) 

where h nj f m and g n _ m are discrete functions and n and m are 
discrete variables. 


1.3-2 CONVOLUTION THEOREM 

Convolution is a comparatively complicated operation because 
integration is involved. By using the Convolution Theorem we are able 
to simplify the calculation of convolution greatly. 

The Convolution Theorem states that the convolution operation in one 
domain corresponds to a complex multiplication operation in the other 
domain (Bracewell, 1984). 

If there is a convolution operation in the function domain as given 
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in equation (1.5) above, then the Fourier Transforms F(s), G(s) and 
H(s) have the following multiplicative relation: 

F(s)G(s) = H(s) . (1.7) 

To use the Convolution Theorem to calculate the convolution between 
f(x) and g(x) : 

h(x) = f(x)*g(x) , (1 .8) 

first we take the transforms of f(x) and g(x), which are F(s) and 
H(s), then rautiply the transforms F(s) and G(s): 

H(s) = F(s)G(s). (1 .9) 

The convolution result is obtained by simply talcing the inverse 
transform of the product of F(s) and G(s). Thus the operation in the 
transform domain can be faster to obtain on the computer. 


1.3.3 DECONVOLUTION 

The relationship between the output signal of a linear 
shift-invariant instrument and its input signal has been established. 
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The input convolved with the impulse response function of uhe 
instrument gives the output signal. But it is often the situation that 
we know the output signal and the impulse response function of the 
instrument. What we need to find is the input signal. In o her words 
we need to recover the actual input signal f in terms of the response 
function g and the output signal h. This process is called 
deconvolution. 

There are many deconvolution methods. One of the most popular 
methods is inverse filtering. The following diagram shows inverse 
filtering or deconvolution in the function domain. 


h 


f 


> > 


g-1 ; > > — 

1 


h = f*g, 


( 1 . 10 ) 


f = h*g~1 = f*g*g-1 , (1.11) 
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By the Convolution Theorem the above operations in the transfom 
domain become multiplication and division: 


F(s)G(s) = H(s) (1.12) 

and 


F(s) = H(s)/G(s). (1.13) 


This result is undefined where G(s) = 0. Tne principal solution 
suggested by Bracewell and Roberts (1954) is: 


H(s)/G(s) [ s: G(s)#0 ] 

F p(s) = 

^0 [s:G(s)=0]. 


(1.14) 


Figure 1.1 and 1.2 show the magnitudes of the transforms. Because of 

symmetry, only values at positive frequencies are presented. Figure 

1.1 shows 'Fj, |G| , and ]Hj , while Figure 1.2 shows Jp i i n i . 

pi> i i ana in i . 

In the transform domain, at regions where G(s) is zero, F f a \ ,- a * 

" p\ S ) lto S0 u 

to zero. This is sometimes equivalent to mutiplying F ( s ) a 
truncating function ( rectangular window ) : 
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Magnitude 


Figure 1.2 


i 
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(1.15) 


F p(s) = F(s)n(s/2s c ), 


where II (s/s ) 

V * 


is "the rectangular window with -the cutoff frequency 


Since sinc^s^) i3 e inve^e Fourier Transform of II(s/2s c ), *jy 
ihe Convolution Theorem the principal solution for this case is: 


p(x) = f(x) * s c sinc 2s c x, (1.16) 

in the function domain. The superposition of these effects 

(convolution with the sine inverse transform) is the approximate 

UuLon fp(x). Proper tapering of the rectangular window reduces the 
oscillations introduced by the sine. 

Notice that f - 

pvxj is no u a unique solution to the deconvolution 
problem. Any da^a naving a transform which is zero at regions where 

G(s) is non-zero and is non-zero at regions where G(s) equals zero can 
be added to f(x) to give a solution to the same deconvolution problem 
(Brae ewe 11 and Roberts, 1954; Ioup and Ioup, 1983). 

If additive noise n(x) is present in the data, the output becomes: 

d(x) = h(x) + n(x). (1.17) 
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The resulting experimental data d(x) are: 

d(x) = f(x)*g(x) + n(x) . (1.18) 

If no noise removal is applied, the principal solution becomes: 

F p(s) = H(s)/G(s) + N(s)/(r(s) (1.19) 

Beyond s Cj -the cutoff frequency of G(s), H(s) is zero. However, the 
noise N(s) is not necessarily zero in this region. For this region, we 
can simply remove the noise entirely first because H(s) is expected to 
be zero in this region. 

1 .4 ITERATIVE DECONVOLUTION METHODS 

For application to real data, iterative deconvolution methods are 
widely used. The ones involved in this research are: 

1. Van Cittert's Method (1931), and 

2. Always-Convergent Deconvolution Method of 
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Ioup (1981), which will be discussed in the next chapter. 


1-4.1 Van Cittert ' 3 Deconvolution Method 

In 1931 Van Cittert created an iterative deconvolution method. The 

iterations in the function domain can be described try the following 
equations : 

f 0(x) = h(x) 

■^l(x) = fo(x) + [ h(x) - fo(x)*g(x) ] 

■^2(x) = fl(x) + [ h(x) - fi (x)*g(x) ] 

• • . 

• • . 

^n-l(x) = fn-2(x) + [h(x) - fn-2(x)*g(x) ] 

^n(x) = fn-1 (x) + [ h(x) - fn-1 (x)*g(x) ] (1.20) 

Van Cittert recognized "hat the image data h(x) could be considered 
as a first approximation, f o(x)> desired f(x). Then he took the 

difference h(x) - f 0 ( x )*g(x) to be the correction for the second 
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iteration. A similar correction is obtained for each iteration. In 
other words, for each iteration, the approximation result f ( x ) ia -^ e 

sum of the previous approximation f n _, ( x ) correction factor 

which is the difference between the image data h(x) and the 
convolution of the previous approximation f n _ 1 ( x ) with the response 
function g(x). 


In the transform domain the iterations become: 


F 0(s) = H(s) 

P l(s) = Ho(s) + [ H(s) - I'o(s)G(s) ] 

P 2(s) = Pi(s) + [ H(s) - F 1 ( s )G(s) ] 

• • • 

• • » 

* • • 

F n-1 (s) = F n _ 2 (s) + [ H(s) - F n _ 2 (s)G(s) ] 

P n(s) = F n _i (s) + [ H(s) - F(s) n _iG(s) ] (1.21) 


The last equation may also be written as: 


F n = [ 1 + ( 1 — G) 1 + (1-G) 2 + + (i_G)n ] H 
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= { [ 1 - (I-G)IH-I ] / G } H 


( 1 . 22 ) 


Tnis result was obtained by successive substitutions. Prom equation 
(1.22), we can see that in order for the iterations be be convergent, 
Hie following conditions have to be satisfied: 


! i-o ! < i 

for regions where G(s) is non-zero, 

(1.23) 

H(s) = 0 

for other regions. 

(1.24) 


These are the Convergence Criteria (Ioup and Ioup, 1983). In the 
absence of noise, the convergence depends heavily on the impulse 
response function g(x). Bracewell and Roberts (1954) gave those 
conditions originally. Hill (1973) and Hill and Ioup (1976) used ^hese 

conditions to find necessary restrictions on the shape and location of 
the impulse response function g(x). 

In satisfying condition (1.23), the origin of the impulse res^nse 
function g(x) is of importance (Hill, 1973; Hill and Ioup, 1976) ! 
Sometimes we can adjust the origin of the impulse response function 
g(x) to such a position that G(s) satisfies condition (1.23). 
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CHAPTER II 


ALWAYS-CONVERGENT DECONVOLUTION METHOD OF IOUP 

Van Cittert's Iterative Deconvolution Method is a very effective 
method and has been used since the early thirties. But in order to use 
this method, the impulse response function has vo satisfy Convergence 
Criteria (1 .23 and 1 .24), which means that the transform of the 
impulse response function has to he inside the uni o circle cen tiered on 
1+iO in the complex plane. Many impulse response functions do not 
satisfy Convergence Criteria (1.23 and 1 .24), and often there is nou 
too much we could do to mal® them satisfy (1 .23 and 1 .24) • In 1981 , 
Ioup created a new method to overcome the convergence problem, which 
is the so-called Always-Convergent Iterative Deconvolution Method 
(Ioup, 1 981 , Whitehorn 1980, Ioup and Whitehorn, 1981, Amini 1987). 

The Always-Convergent Method of Ioup is a modified version of van 
Cittert's method. Two modifications are made. First, the transform of 
the imp uls e response function G(s) is replaced by 


^m(s) = !G(s)i / lG(s)imax> 


( 2 . 1 ) 
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and secondly, in "Che iteration equations, H(s) is replaced by the 
product of ihe principal solution and G m ( s ). 

Prom t he definition of G m (g) f we can see that G m (s) is always real 
and between 0 and 1 , which guarantees that G m ( a ) satisfies Convergence 
Criterion (1.23), thus the iterations are always convergent. 

Applying the above modifications to the van Cittert iterations in 
the transform do m ain, the Always-Convergent Method in the transform 
domain can be described by the following equations: 


p 0= H 

F 1 - F 0 + [ Fp - P 0 ] S B 

• • • 

• • • 

• • • 

F n = F n-1 + [ F p ~ F n-1 1 Gtn. (2.2) 

With successive substitutions, the above equations can be written in 
the general form: 


F n = [ 1 - (1-G)(1-G m ) n ] H/G. (2.3) 
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For regions where G(s) is zero, then: 



( 2 . 4 ) 


The Always-Convergent Deconvolution Method can be also used in the 
time domain. The iterations are given by the following equations: 


f 0 = ^0 

1 1 = To + C f p - f 0 ] * g m 

• • • 

• • • 

• • « 

f n * f n-1 + [ f p ~ T n -1 ] * gm, (2.5) 


w h ere Sqj and fp are the inverse Fourier Transforms of G m and Fp. 

With the Always-Convergent Method in the time domain, since the 
convolution operation is involved, it is harder to do the iterations 
and takes more computer time, but one usually has more information 
about the original signal than the transform. Thus during the 
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deconvolution process, one may be able to check toe result and make 
corrections to toe result to make toe deconvolution more precise. For 
example, for a positive signal, one knows that toe result should be 
non-nega „ive, thus toe non-negative constraint may be applied to the 
deconvolved signal to improve toe deconvolution. 


For toe Always-Convergent Method in toe transform domain, only- 
complex multiplication is involved, so it is easier and much faster to 
deconvolve on toe computer. The disadvantage is that one usually has 
less information about toe transform toan about toe original signal, 
tous it is harder to make adjustments during toe deconvolution process 
to get a more accurate deconvolved result. 

Au regions where G(s) is zero, toe deconvolution equation (2.3) 
becomes : 


[P n(s)][0] = H(s) . 


( 2 . 6 ) 


Obviously any arbitrary function will satisfy (2.6), and we are not 
be able to recover toe input. For all of toe above deconvolution 
metoods, toe solution F(s) is defined to be zero in tois region. This 
does not really solve toe problem, as we still have no idea what toe 
real input is. Howerer, toe application of function domain constraints 
can extend toe solution into tois region ! 
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CHAPTER III 


NOISE 


In the real world there will be always noise present. Data 
processing systems are usually required to handle a large assortment 
of signals in the presence of noise. There are many iypes of noise: 
noise could be ihe wrong signal, or the right signal out of place, and 
what is noise to one observer might be signal to another. In any 
application, care must be given to the classification of which signals 
are considered to be useful and which ones are considered undesirable, 
or noise (Ronbinson, 1980). 

The noise studied in this research is Gaussian distributed noise. 
Gaussian noise has a Gaussian or bell curve distribution about any 
given data point. If enough noise cases are generated, a plot of 
occurrence versus the magnitude of that amplitude approaches a 
Gaussian shape (Pig.3.1 ). 
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There are two types of Gaussian distributed noise: constant and 
ordinate-dependent noise. For constant Gaussian noise, the width of 
the bell curve is fixed; it is Gaussian noise with the same standard 
deviation at each point. However, ordinate-dependent Gaussian noise is 
that for which the width of the bell curve depends on the ordinate 
value of the data point. 

To add Gaussian noise to the data, one can use the well-known result 
that the sum of comparatively few random numbers from a uniform 
distribution gives a very good approximation to a normal Gaussian 
distribution (H ammi ng, 1962). Utilization of the Central-Limit Theorem 
leads to the above conclusion (Hamming, 1962). On the computer, in the 
case of a decimal machine (expressed as a base 10 number) it is 
customary to use 12 random numbers to get a variance of 1 (Hamming, 
1962). Since the mean in this case is not zero (the mean is half of 
the number of points), a necessary amount (in this case 6) must be 
subtracted from the sum of the numbers to makp the mean zero. For the 
constant noise case the expression is: 


h ncd) = (Ai - 6) SF 1 / 2 + h(l) (3.1) 

/ 

where the are the sum of random numbers, SF is the scale factor, 
and h is the noise free data. With the help of the SF one can generate 
a noisy data set which has approximately the signal to noise ratio 
(SNR) of interest (Leclere, 1984). 
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CHAPTER IV 


SINGLE-FILTER APPLICATION OF THE ALWAYS-CGNVERGENT 


DECONVOLUTION METHOD TO PHYSICAL DATA 


In uhis chapter, the use of the Always-Convergent Deconvolution 
Method to deconvolve synthetic data is discussed. The input used here 
consists of three narrow Gaussian peaks which are shown in Fig. 4.1. 
The second and third peaks are located close enough together so that 
after convolution there is some overlap. To get the output data, the 
input is convolved with both narrow and wide Gaussian impulse response 
functions. The narrow Gaussian impulse response function contains 9 
poinus as shown in Fig. 4. 2, whale the wide Gaussian contains 21 points 
as shown in Fig. 4 . 3- Gaussian- type noise with various signal- to-noise 
ratios is added wo the output data. The deconvolved result is compared 
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to the input, and the mean square error (MSE) between the deconvolved 
result and "hie input is calculated as a quality measure. The 
minimization of the mean square error between the deconvolved result 
and the known synthetic input is used to determine the optimum 
iteration number. 

The mean square error of the deconvolved result F' and the input F 
in the transform domain is given by: 


MSE = [ 1 ( F'i-Pi )2 ]l/2, 

i= t 


* 


(4.1) 


which is the square root of the sum over all frequencies of the 
squares of the difference between the transform of the deconvolved 
result and the transform of the input. 

Two Gaussian impulse response functions were used since many 
instruments have Gaussian impulse response functions. We use both 
narrow and wide Gaussian impulse response functions, expecting to 
cover approximately the range of typical instruments. 

Discrete Fourier Transforms of the sampled impulse response function 
g(x) and the output h(x) (with or without noise) are taken for both 
narrow and wide Gaussian impulse response functions, using an FFT3D 
subroutine (DEL). Equation (2.3) is used to perform the deconvolution 
iterations in the transform domain. 
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The mean square error between the known synthetic input and the 
deconvolved result is calculated for each iteration and compared to 
that of the previous iteration. Previous research results (Amini, 
1986) show that the curve of mean square error versus iteration number 
under certain signal- to-noise ratios has the shape of Fig. 4. 7. This 
means that the mean square error decreases as the iteration number 
increases until the optimum iteration is reached. Beyond the optimum 
iteration number, the mean square error increases as the iteration 
number increases. During the process of deconvolution, once the MSE is 
found to increase compared ~so the previous iteration, we simply 
terminate the deconvolution and set that iteration number to be the 
optimum iteration number. 

Table 4.1 shows the optimum iteration number for the narrow Gaussian 
case at a signal to noise ratio of 10 to 150. The optimum iteration 
number is 3 at a SNR of 10, and increases to about 19 at a SNR of 150; 
the optimum iteration number is larger at lower SNR regions and 
smaller at higher SNR regions. 

Table 4.2 shows similar values for the wide Gaussian case. The 
optimum iteration number increases from 3 at a SNR of 10 to 53 at a 
SNR of 150. This has the same behavior as that of the narrow Gaussian 
case, except the magnitudes are larger. 

At the regions of low signal to noise ratio, the noise is large, and 
the large noise will terminate the iterations quickly. That is why the 
optimum iteration numbers are smaller in this region. As the signal to 
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noise ratio increases, the noise becomes smaller and more signal can 

* rea ^ red ' U Wl11 iterations to reach *e optimum level 

of restoration, thus the optimum iteration number becomes larger in 
"this region. 


If we compare the optimum iteration numbers for these two cases, we 
eee that the narrow Gaussian case has a lower value than the wide 

Gaussian. This is a reasonable result, because the narrow Gaussian 

impulse response function has a wider transform than that of the wide 

Gaussian. Thus at most frequencies the G m of . - . 

m oi a narrow Gaussian is 

closer 10 one (or il-G i •, ^ 

mi 13 closer oO 0) than that for the wide 

Gaussian. Prom e,uatio„(2.3) , the closer |M.j l3 ^ ^ ^ 

“ “' Ker « enoe ^ ad Thus the narrow Gaussian reaches the optima 
iteration faster, or has a smaller optimum iteration number. 

Table 4.3 shows the mean square error for the narrow Gaussian case 

“ frequencjr domain at a SSR range from to to 150. To save 
computer time, the mean square error is calculated in the transform 
domain. The mean square error decreases rapidly from 247.7 at a SNR of 

*° 34 ' 9 a “ a ^ of ,5 °- The data show clearly that the mean square 
error is ve^r large at regions of low signal to noise ratio and drops 
rapidly as the signal to noise ratio increases. This is obviously a 
reasonable result. Under the same input signal, the mean square error 
depends on two factors: noise and the impulse response function. At 
rations of low si^ to noise ratio, the noise is large, and we 
simply get a large mean square error due to large noise. At regions of 
high signal to noise ratio, the noise is small and the impulse 
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response function becomes a more important source of mean square error 
due to lack of restoration. 

The mean square error depends heavily on the shape of the imp uls e 
response function. Figure 4-4 shows the transform of a narrow Gaussian 
impulse re&ponjae func ^ion . We see that the transform has only one 
almost zero point and increases rapidly to some considerably large 
amplitudes on both sides. In other words, at almost all frequencies 
the transform G(s) of the impulse response function has relatively 
large amplitudes. If G(s) is non zero, we do not have to assume H(s) 

be zero. If the amplitude of G(s) is not very small, a little noise 
added to H(s) will not affect the result of H(s)/G(s) significantly. 
Thus the mean square error is very small at high signal- to-noise 
ratios for the narrow Gaussian case. 

Table 4.4 presents similar MSE results in the transform domain for 
the wide Gaussian case. The mean square error decreases slowly from 
559*0 ao a SNR of 10 to 449*2 at a SNR of 150. The mean square error 
is very large in the low SNR region (high noise region) and is only 
slightly smaller in the high SNR region (low noise region). Tnis is 
obviously a different situation from the narrow Gaussian case. At 
regions of low signal to noise ratio (high noise region) , it is 
natural that due to "the large noise, the mean square error is large. 
But in the high SNR region (low noise region), the mean square error 
is still very large, only slightly smaller than that of the high noise 
region. Consider the impulse response function. Figure 4.5 is the 
transform of the wide Gaussian impulse response function. In a very 
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wide frequency region, the transform appears to be zero. Actually the 
transform is not zero but is very small (on the order of 10~7 or 
10" 6 )* Thus ^en we do the division H(s)/G(s), a small amount of noise 
will cause significantly large mean square errors. 

From the above results, we see that the narrow Gaussian case has a 
smaller mean square error and thus a better deconvolution result; the 
wide Gaussian has a larger error and worse deconvolution results. 

Wraparound error can also be a significant error source during the 
deconvolution process. If the inverse of the impulse response function 
and the output signal are not narrow enough compared to the whole time 
interval, errors will arise at the edges due to overlapping effects. 
This is called the wraparound error. One way to reduce or get rid of 
the wraparound error is to add more zero points to the original 
function to increase the length. But the question is: will this affect 
the deconvolution result ? BenSueid (1987) showed that when no noise 
is present, adding zeros to he original function will not affect he 
deconvolution result if he means square error is calculated in he 
time domain. later we will show hat even when noise is present, we 
will have he same result in he time domain. Thi 3 means hat in order 
to reduce he wraparound error one can add as many zeros as he wants 
t» he original function. But if zeros are added, it will take n»re 
computer time. So one needs to choose an optimum lengh (adding he 
proper number of zeros) which will minimize he wraparound error but 
do slow down he computer unnecessarily. 
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Table 4.5 shows the optimum iteration numbers and mean square errors 
(is tee time domain) for different lengths for tee narrow Gaussian 
case at a SNR of 10 to 150. The optimun iteration numbers and mean 
square errors change significantly when tee length increases from 52 
points to 64 points; after teat, they do not change at all while tee 
length keeps increasing from 64 to 128 and from 128 to 256 and so on. 
Prom tee fact teat tee optimum iteration numbers and mean square 
errors do not change when tee length takes a value longer than 64 
points, we can conclude teat adding seroe beyond 64 points to tee 
original taction even when noise is present does not affect tee 

deconvolution result. Ms conclusion allows one to reduce tee 
wraparound error for the noise situation. 


The mean square errors and tee optimum iteration numbers change 
significantly when tee length increases from 52 points to 64 points 
and remain constant when tee tatfh is »er 64 points for teis case.’ 
Therefore, tee wraparound error is significant when tee length is 52 
points and vanishes when tee length is over 64 points. Thus we can 
Simply Choose 64 to be tee optimum length for teis case to reduce 
wraparound error and save computer time. 


Table 4.6 shows similar quantities for tee wide Gaussian case. The 
results have tee same tehavior as tee narrow Gaussian case. Thus tee 
same conclusion can be made as for tee narrow Gaussian case. 

One important thing the above results show is teat tee narrow 
Gaussian case has better deconvolution results than tee wide Gaussian 
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case. But does a narrow impulse response function always have a better 
deconvolution result than a wide impulse response ? If so, why ? What 
can we do to get a better deconvolution result ? Analysis in "the 
transform domain has been made on whether a deconvolution is 
successful ! How about analysis in the time domain ? Curiosity mates 
me look for the answers to these questions ! 


Figure 4.12 illustrates some general impulse response functions 
(sometimes called windows). Curve g c is an given impulse response 
function. As curve g c becomes wider yet, it takes the shape of curve 

®b* As it gets wider, it becomes a flat "curve" (curve ga), which is a 
norizon ual straight line. On the other hand, if curve g gets 
narrower, it becomes a delta function (curve g d ). Therefore this graph 
shows windows of a range of widths, from the narrowest delta fine t ion 
(curve g d ) -fo -th e widest horizontal line (curve g a ). 


Figure 4.13 is an input signal f with two peaks. To see what happens 
in the time domain during the deconvolution process, this input f is 
convolved with each of the above windows individually. In figure 4.14, 

n a» hb, be 3/1(1 h d are convolutions of f with ga, g^, g c and g d . 


Recall that the convolution can be defined by the following: 


h(x) = 



f(u) g(x-u) 


du 


(4.2) 


To get the convolution value at a certain point x, the window is 
reversed and centered on x, then multiplied by the input and 
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integrated ewer the period (to avoid wraparound errors, the period is 
selected much larger than the width of the windows). This is actually 
an average", calculation of the input with a "weight”, function g(x-u). 
If the window g(x) is narrow, then we are "averaging", the input in a 
narrow region for each convolution value. Thus we will loose less 
detail about the input during the convolution process and the 
convolution will have a shape more like the input. The deconvolution 
will be easier in this situation. On the other hand, if we have a wide 
window, then we are "averaging”, the input over a wide region for each 
convolution value. More details of the input get lost during 

convolution and it will be harder to recover the signal. This can be 
observed in Figure 4.14- 


In Figure 4.12, g d l3 a delta function or the narrowest window, hd 
is the convolution of g, witt tte input, which is the same as tte 
input because the delta function is the identity for convolution. When 
the input is convolved with a delta function, it does not change at 

all. It is easy to do such a deconvolution when the output is 
identical to the input. 


Now consider the windows g ^ . 

s c a™ «b. Window g c 13 narrower than 

window g b . The outputs ^ and h b are the convolution of g c and g b with 

the input f. We see that h Q lies between hd (the same as the i n p u t) 

' n bt which means hr» has values rO^vao-n -s_ - , 

c v ^ues closer uo 1 he input or carries more 

information about the inpit ttan h b . Thus it is easier * ^ 

inpio from h Q (the result of convolution with the narrower window) 
tnan from h b (the result of convolution with the wider window). 
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If the window keeps getting wider, finally it will become the 
horizontal line g a . The output h a ( convolution of ga with the input 
f) is a constant; indeed, the convolution of this window with any kind 
of inpuu will be a constant. This is not a surprising result, because 
we are actually "averaging", the entire input. Any input convolved with 
this window becomes a constant output. Obviously the solution to this 
deconvolution problem can be any function, or in other words, we are 
not able to get the solution we wanted. 

We can also make a similar analysis in the transform domain and 
reach the same conclusion. 

Figure 4-. 15 illustrates the transforms of the windows shown in 
Figure 4.12. G &> g c and G d are the transforms of g a , gc and 

S(j. Note that the transform of a delta function is a constant and the 
transform of a constant is a delta function. Also we see that the 

narrower window h c has a wider transform than that of the wider window 

V 


For the convolution with a delta function, deconvolving in the 
transform domain is simply a* division H^. Slnce 0(J l3 a 00Mltantf 

this is almost as simple as the parallel deconvolution in the time 
domain. 

Comparing the transforms G ru „„ ___ j. A . 

c we see that the transform Gq of 

the narrower window (g c ) ia wider than the transform of the wider 
window (g b ) . This means G^ has more zero points or more small 
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magnitude points uhan G^ m t 0 <j 0 deconvolution, one has to perform 
Ihe division H/G. If in a certain region G is zero, then the division 
is undefined. In other words, the information in the input in this 
region gets lost. If G is not zero but has very small magnitude, then 
as discussed before, large errors will be produced in this region. 
Therefore from the transform domain the narrower window will retain 
more details about the input and produce less error. 

The transform of the constant function g a ia G a , which is a delta 

function. H a is required to be zero in -the entire transform domain 
except for one point. In other words, all the information about the 
input nas been lost in the transform domain except the value at the 
origin. Any input F will give a solution to the deconvolution problem. 
This the same conclusion as was reached in the time domain. 
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CONCLUSION 


In "this thesis, the Always-Convergent Iterative Noise Removal and 
Deconvolution Method of Ioup (Ioup, 1 981 ) is applied as a single filter 
oO deconvolution with both narrow and wide Gaussian impulse response 
functions. The wraparound error for both cases is also studied. The 
results show that the narrow Gaussian impulse response converges 
faster and has smaller mean square error than the wide Gaussian 
impulse response. In general, deconvolution using a narrow window 
usually produces better results (converging faster and having smaller 
mean square error) than a wide window. In our particular case, "the 
deconvolution result of the narrow Gaussian case matches the original 

input very well, while the result using the wide Gaussian does not 
match the input so well. 

In our case, the narrow Gaussian impulse response has smaller 
optimum iteration numbers and smaller mean square error than the wide 
Gaussian, which matches the above conclusion. 

This thesis also shows that for a deconvolution problem, even when 
noise is present one can still add zeros to the original function 
(increasing the length) to reduce the wraparound error. This result 
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enables one to minimize toe wraparound error when noise is present. In 
order to save computer time, a optimum length should be chosen. To 
obtain toe optimum length, one has to minimize both toe wraparound 
error and toe computer time. In toe narrow Gaussian case, toe optimum 
iteration numbers and toe mean square errors change significantly when 
tte lengto increases from toe original 32 points to 64 points and 
remain unchanged when toe lengto is over 64 points. This means toe 
wraparound error is only significant at toe lengto of 32 points and 
vanishes after toe lengto is longer toan 64 points. Thus toe optimum 
length is obviously 64 points for this case. In toe same way, toe 
optimum lengto for toe wide Gaussian case is chosen to be 88 points. 
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TABLE 4.2 


OPTIMUM ITERATION NUMBERS OF THE WIDE GAUSSIAN CASE FROM 


SIGNAL TO NOISE RATIO OF 10 TO 150 

( 50 CASES ) 


SIGNAL TO NOISE RATIO 

OPTIMUM ITERATIONS 

10 

3.0 

20 

6. 1 

30 

9.8 

40 

12.9 

50 

17.0 

60 

20.8 

70 

23.9 

80 

27.7 

90 

31.8 

100 

35.2 

1 10 

39.0 

120 

43. 1 

130 

46.2 

140 

50. 1 

150 

53.0 
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TABLE 4.1 


OPTIMUM 


TERAT ION NUMBERS OF THE NARROW GAUSSIAN CASE FROM 


SIGNAL TO NOISE RATIO OF 10 TO 150 

( 50 CASES ) 


SIGNAL TO NOISE RATIO 
10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
110 
120 
130 
140 
150 


OPTIMUM ITERATIONS 

2.8 

4.0 

5.9 

7.8 

9.9 

11.0 
12.2 

13.1 

14.0 

15.1 
16.0 

16.9 
18.0 

18.9 

19.2 
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TABLE 4.4 


THE MEAN SQUARE ERROR OF THE WIDE GAUSSIAN CASE FROM 


SIGNAL TO NOISE RATIO OF 10 TO 150 
( 50 CASES ) 


SIGNAL TO NOISE RATIO 

MEAN SQUARE ERROR 

10 

559.0 

20 

531.1 

30 

515. 1 

40 

502.9 

50 

493.4 

60 

485.5 

70 

478.8 

80 

473.2 

90 

468.3 

100 

464. 1 

1 10 

460.4 

120 

457.1 

130 

454.2 

140 

451 .5 

150 

449.2 
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TABLE 4.3 


THE MEAN SQUARE ERROR OF THE NARRCW GAUSSIAN CASE FROM 

SIGNAL TO NOISE RATIO OF 10 TO 150 

( 50 CASES ) 


SIGNAL TO NOISE RATIO 

MEAN bUUMnc tnn 
247.7 

10 

159.5 

20 

122 . 5 

30 

100.9 

40 

86.3 

50 

75.3 

60 

66.8 

70 

60.0 

80 

54.5 

90 

49.8 

100 

45.9 

110 

42.6 

120 

39.7 

130 

37.2 

140 

34.9 

150 
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TABLE 4.5 


MSE RESULTS AT DIFFERENT LENGTHS, NARROW GAUSS I AN 


SNR 


LENGTH 



1 32 

1 64 i 

1 128 

1 256 1 

1 


MSE 


1 

10 

82 

.7 

80.3 

80.3 

80.3 

20 

56 

.6 

53.6 

53.6 

53.6 

30 

41 

.7 

38.8 

38.8 

38.8 

40 

32. 

.3 

30.1 

30. 1 

30. 1 

50 

26, 

.5 

24.7 

24.7 

24.7 

60 

22 , 

.4 

21.0 

21 .0 

21 .0 

70 

19. 

.4 

18.2 

18.2 

18.2 

80 

17. 

1 

16.0 

16.0 

16.0 

90 

15. 

2 

14.3 

14.3 

14.3 

100 

13. 

7 

12.9 

12.9 

12.9 

1 10 

12. 

5 

11.8 

11 .8 

11.8 

120 

1 1 . 

5 

10.8 

10.8 

10.8 

130 

10. 

6 

10.0 

10.0 

10.0 

140 

9. 

9 

9.3 

9.3 

9.3 

150 

9. 

2 

8.7 

8.7 

8.7 
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